function [jobs, jobhelps, funcfiles] = spm5_preprojobs(sf, fp, ap, opts)
% spm5_preprojobs  - create jobs and jobhelps from options
%
% FORMAT:       [jobs, jobhelps, ff] = spm5_preprojobs(sf, fp, ap, opts)
%
% Input fields:
%
%       sf          subject (base) folder
%       fp          functional folder pattern (e.g. f*/r*)
%       ap          anatomical folder pattern (e.g. a*/*)
%       opts        1x1 struct with options
%        .afp       anatomical file pattern (default: '*.img' / '*.nii');
%        .ffp       anatomical file pattern (default: '*.img' / '*.nii');
%        .jobfile   save output into a mat file (full or partial filename)
%        .jobrun    flag, run jobs?
%        .rinterpe  realignment estimation spline degree (default: 2)
%        .rinterpr  realignment reslicing spline degree (default: 4)
%        .rqual     realignment quality (accuracy trade-off default: 0.9)
%        .rrtm      realignment register-to-mean flag (default: true)
%        .rsep      realignment sampling grid resolution (default: 3mm)
%        .smk       smoothing kernel (either 1x1 or 1x3 double, default: 6)
%        .st        flag, perform slice time correction (default: true)
%        .sto       slice acquisition order, either of
%                   'a', 'aie', {'aio'}, 'd', 'die', 'dio' or numeric list
%        .str       reference slice (default: 1)
%        .stta      TA (default: auto)
%        .sttr      TR (default: 2 sec)
%        .wbb       normalize writing bounding box
%                   (default: [-78, -114, -54; 78, 78, 84])
%        .winterp   normalize writing interpolation spline (default: 1)
%        .wvox      normalize writing voxel size (default: [3, 3, 3])
%
% Output fields:
%
%       jobs        cell array with jobs suitable for spm_jobman
%       jobhelps    cell array with helps (required when saved as jobs.mat)
%       ff          cell array with functional files (after preprocessing)

% Version:  v0.8a
% Build:    9102621
% Date:     Oct-26 2009, 9:18 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
   ~ischar(sf) || ...
    isempty(sf) || ...
    exist(sf(:)', 'dir') ~= 7 || ...
   ~ischar(fp) || ...
    isempty(fp) || ...
   ~ischar(ap) || ...
    isempty(ap)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid or missing argument supplied.' ...
    );
end
sf = sf(:)';
fp = fp(:)';
ap = ap(:)';

% check options
if nargin < 4 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'afp') || ...
   ~ischar(opts.afp) || ...
    isempty(opts.afp)
    opts.afp = {'*.img', '*.nii'};
else
    opts.afp = opts.afp(:)';
end
if ~isfield(opts, 'ffp') || ...
   ~ischar(opts.ffp) || ...
    isempty(opts.ffp)
    opts.ffp = {'*.img', '*.nii'};
else
    opts.ffp = opts.ffp(:)';
end
if ~isfield(opts, 'jobfile') || ...
    isempty(opts.jobfile) || ...
   ~ischar(opts.jobfile)
    opts.jobfile = '';
end
if ~isfield(opts, 'jobrun') || ...
    numel(opts.jobrun) ~= 1 || ...
   (~isa(opts.jobrun, 'double') && ...
    ~islogical(opts.jobrun))
    opts.jobrun = false;
end
if ~isfield(opts, 'rinterpe') || ...
   ~isa(opts.rinterpe, 'double') || ...
    numel(opts.rinterpe) ~= 1 || ...
    isinf(opts.rinterpe) || ...
    isnan(opts.rinterpe) || ...
   ~any((1:7) == opts.rinterpe)
    opts.rinterpe = 2;
end
if ~isfield(opts, 'rinterpr') || ...
   ~isa(opts.rinterpr, 'double') || ...
    numel(opts.rinterpr) ~= 1 || ...
    isinf(opts.rinterpr) || ...
    isnan(opts.rinterpr) || ...
   ~any((1:7) == opts.rinterpr)
    opts.rinterpr = 4;
end
if ~isfield(opts, 'rqual') || ...
   ~isa(opts.rqual, 'double') || ...
    numel(opts.rqual) ~= 1 || ...
    isinf(opts.rqual) || ...
    isnan(opts.rqual) || ...
   	opts.rqual <= 0 || ...
    opts.rqual > 1
    opts.rqual = 0.9;
end
if ~isfield(opts, 'rrtm') || ...
    numel(opts.rrtm) ~= 1 || ...
   (~isa(opts.rrtm, 'double') && ...
    ~islogical(opts.rrtm))
    opts.rrtm = true;
elseif isa(opts.rrtm, 'double')
    opts.rrtm = (opts.rrtm ~= 0);
end
opts.rrtm = 0 + opts.rrtm;
if ~isfield(opts, 'rsep') || ...
   ~isa(opts.rsep, 'double') || ...
    numel(opts.rsep) ~= 1 || ...
    isinf(opts.rsep) || ...
    isnan(opts.rsep) || ...
   	opts.rsep < 1 || ...
    opts.rsep > 8
    opts.rsep = 3;
end
if ~isfield(opts, 'smk') || ...
   ~isa(opts.smk, 'double') || ...
   ~any([1, 3] == numel(opts.smk)) || ...
    any(isinf(opts.smk) | isnan(opts.smk) | opts.smk < 0 | opts.smk > 16)
    opts.smk = 6;
end
if numel(opts.smk) == 1
    opts.smk = opts.smk([1, 1, 1]);
end
if ~isfield(opts, 'st') || ...
    numel(opts.st) ~= 1 || ...
   (~isa(opts.st, 'double') && ...
    ~islogical(opts.st))
    opts.st = true;
end
if isa(opts.st, 'double')
    opts.st = (opts.st ~= 0);
end
if ~isfield(opts, 'sto') || ...
   (~ischar(opts.sto) && ...
    ~isa(opts.sto, 'double'))
    opts.sto = 'aio';
end
if ischar(opts.sto)
    if ~any(strcmpi(opts.sto(:)', {'a', 'aie', 'aio', 'd', 'die', 'dio'}))
        opts.sto = 'aio';
    else
        opts.sto = lower(opts.sto(:)');
    end
else
    opts.sto = opts.sto(:)';
    if any(isinf(opts.sto) | isnan(opts.sto) | opts.sto < 1 | opts.sto ~= fix(opst.sto)) || ...
        numel(opts.sto) ~= numel(unique(opts.sto))
        warning( ...
            'BVQXtools:BadArgument', ...
            'Invalid slice order argument given.' ...
        );
        opts.sto = 'aio';
    end
end
if ~isfield(opts, 'str') || ...
   ~isa(opts.str, 'double') || ...
    numel(opts.str) ~= 1 || ...
    isinf(opts.str) || ...
    isnan(opts.str) || ...
    opts.str < 1 || ...
    opts.str ~= fix(opts.str) || ...
   (isa(opts.sto, 'double') && ...
    opts.str > max(opts.sto))
    opts.str = 1;
end
if ~isfield(opts, 'stta') || ...
   ~isa(opts.stta, 'double') || ...
    numel(opts.stta) ~= 1 || ...
    isinf(opts.stta) || ...
    isnan(opts.stta) || ...
    opts.stta <= 0
    opts.stta = [];
end
if ~isfield(opts, 'sttr') || ...
   ~isa(opts.sttr, 'double') || ...
    numel(opts.sttr) ~= 1 || ...
    isinf(opts.sttr) || ...
    isnan(opts.sttr) || ...
    opts.sttr <= 0
    opts.sttr = 2;
end
if ~isempty(opts.stta) && ...
    opts.stta >= opts.sttr
    opts.stta = [];
end
if ~isfield(opts, 'wbb') || ...
   ~isa(opts.wbb, 'double') || ...
   ~isequal(size(opts.wbb), [2, 3]) || ...
    any(isinf(opts.wbb(:)) | isnan(opts.wbb(:))) || ...
    any(diff(opts.wbb) < 0)
    opts.wbb = [-78, -114, -54; 78, 78, 84];
end
if ~isfield(opts, 'winterp') || ...
   ~isa(opts.winterp, 'double') || ...
    numel(opts.winterp) ~= 1 || ...
   ~any((1:4) == opts.winterp)
    opts.winterp = 1;
end
if ~isfield(opts, 'wvox') || ...
   ~isa(opts.wvox, 'double') || ...
   ~any([1, 3] == numel(opts.wvox)) || ...
    any(isinf(opts.wvox) | isnan(opts.wvox) | opts.wvox < 1 | opts.wvox > 16)
    opts.wvox = 3;
end
if numel(opts.wvox) == 1
    opts.wvox = opts.wvox([1, 1, 1]);
end

% load job specification file
try 
    load([bvqxtools_path('contrib') filesep 'spm' filesep 'spm5_prepro.mat']);
catch
    error( ...
        'BVQXtools:BadFile', ...
        'SPM5 jobs template file not found or corrupt.' ...
    );
end

% set required SPM5 files
try
    spmp = fileparts(which('spm'));
    tpm = { ...
        [spmp filesep 'tpm' filesep 'grey.nii']; ...
        [spmp filesep 'tpm' filesep 'white.nii']; ...
        [spmp filesep 'tpm' filesep 'csf.nii']};
    if exist(tpm{1}, 'file') ~= 2 || ...
        exist(tpm{2}, 'file') ~= 2 || ...
        exist(tpm{3}, 'file') ~= 2
        error('FILE_NOT_FOUND');
    end
    jobs{4}.spatial{2}.preproc.opts.tpm = tpm;
catch
    error( ...
        'BVQXtools:FileNotFound', ...
        'At least one required file of SPM5 was not found.' ...
    );
end

% make settings (without further checks)
jobs{1}.temporal{1}.st.tr = opts.sttr;
jobs{1}.temporal{1}.st.refslice = opts.str;
jobs{2}.spatial{1}.realign{1}.estwrite.eoptions.interp = opts.rinterpe;
jobs{2}.spatial{1}.realign{1}.estwrite.eoptions.quality = opts.rqual;
jobs{2}.spatial{1}.realign{1}.estwrite.eoptions.rtm = opts.rrtm;
jobs{2}.spatial{1}.realign{1}.estwrite.eoptions.sep = opts.rsep;
jobs{2}.spatial{1}.realign{1}.estwrite.roptions.interp = opts.rinterpr;
jobs{3}.util{1}.imcalc.options.interp = opts.rinterpr;
jobs{4}.spatial{3}.normalise{1}.write.roptions.bb = opts.wbb;
jobs{4}.spatial{3}.normalise{1}.write.roptions.interp = opts.winterp;
jobs{4}.spatial{3}.normalise{1}.write.roptions.vox = opts.wvox;
jobs{5}.spatial{1}.smooth.fwhm = opts.smk;

% find folders and files
if any(ap == filesep)
    [afp, ap, ape] = fileparts(ap);
    ap = [ap, ape];
    afp = [filesep afp];
else
    afp = '';
end
anatfolders = findfiles([sf afp], ap, 'dirs', 'depth=1');
if isempty(anatfolders)
    error( ...
        'BVQXtools:FolderNotFound', ...
        'No anatomical folders found...' ...
    );
end
if any(fp == filesep)
    [ffp, fp, fpe] = fileparts(fp);
    fp = [fp, fpe];
    ffp = [filesep ffp];
else
    ffp = '';
end
funcfolders = findfiles([sf ffp], fp, 'dirs', 'depth=1');

% get anatomicals
aff = numel(anatfolders);
if aff == 0
    error( ...
        'BVQXtools:NoFoldersFound', ...
        'No anatomical folders found.' ...
    );
end
for sc = 1:aff
    anatfiles = findfiles(anatfolders{sc}, opts.afp, 'depth=1');
    if ~isempty(anatfiles)
        anatfolders = anatfolders(sc);
        break;
    end
end
if isempty(anatfiles)
    error( ...
        'BVQXtools:NoFilesFound', ...
        'No anatomical file(s) found.' ...
    );
end
anatfiles{1} = [anatfiles{1} ',1'];

% get number of scans and get files in folders
nff = numel(funcfolders);
if nff == 0
    warning( ...
        'BVQXtools:BadArgument', ...
        'No functional folders found.' ...
    );
    funcfiles = cell(1, 0);
    if ~opts.st
        jobs(1) = [];
    end
    return;
end
funcfiles = cell(1, nff);
for sc = nff:-1:1
    funcfiles{sc} = findfiles(funcfolders{sc}, opts.ffp, 'depth=1');
    if isempty(funcfiles{sc})
        warning( ...
            'BVQXtools:NoFilesFound', ...
            'No functional files found in ''%s''.', ...
            funcfolders{sc} ...
        );
        funcfiles(sc) = [];
    end
    if numel(funcfiles{sc}) == 1 && ...
        strcmpi(funcfiles{sc}{1}(end-3:end), '.nii')
        try
            niit = [];
            niit = BVQXfile(funcfiles{sc}{1});
            if ~isBVQXfile(niit, 'hdr')
                error('NO_NIFTI');
            end
            nvol = size(niit.VoxelData, 4);
        catch
            clearbvqxobjects({niit});
            warning( ...
                'BVQXtools:BadFile', ...
                'Invalid NIftI file ''%s''.', ...
                funcfiles{sc}{1} ...
            );
            funcfiles(sc) = [];
            continue;
        end
        niit.ClearObject;
        funcfiles{sc} = funcfiles{sc}(ones(nvol, 1));
        for ic = 1:nvol
            funcfiles{sc}{ic} = sprintf('%s,%d', funcfiles{sc}{ic}, ic);
        end
    else
        for ic = 1:numel(funcfiles{sc})
            funcfiles{sc}{ic} = [funcfiles{sc}{ic} ',1'];
        end
    end
end

% get number of slices and possibly auto settings for ST
tvol = spm_vol(funcfiles{1}{1});
nslc = tvol.dim(3);
if isempty(opts.stta)
    opts.stta = opts.sttr * (nslc - 1) / nslc;
end
if ischar(opts.sto)
    switch (opts.sto)
        case {'a'}
            opts.sto = 1:nslc;
        case {'aie'}
            opts.sto = [2:2:nslc, 1:2:nslc];
        case {'aio'}
            opts.sto = [1:2:nslc, 2:2:nslc];
        case {'d'}
            opts.sto = nslc:-1:1;
        case {'die'}
            fslc = 2 * floor(0.5 * nslc);
            sslc = 2 * floor(0.5 * (nslc - 1)) + 1;
            opts.sto = [fslc:-2:1, sslc:-2:1];
        case {'dio'}
            fslc = 2 * floor(0.5 * (nslc - 1)) + 1;
            sslc = 2 * floor(0.5 * nslc);
            opts.sto = [fslc:-2:1, sslc:-2:1];
    end
end
jobs{1}.temporal{1}.st.nslices = nslc;
jobs{1}.temporal{1}.st.so = opts.sto;
jobs{1}.temporal{1}.st.ta = opts.stta;

% begin to build scans/data into jobs -> slice timing
jobs{1}.temporal{1}.st.scans = funcfiles;
if opts.st
    funcfiles = fnpp(funcfiles, 'a');
end

% realignment
jobs{2}.spatial{1}.realign{1}.estwrite.data = funcfiles;
funcfiles = fnpp(funcfiles, 'r');

% mean image (imcalc)
jobs{3}.util{1}.imcalc.input = cat(1, funcfiles{:});
jobs{3}.util{1}.imcalc.outdir{1} = anatfolders{1};

% coregistration
jobs{4}.spatial{1}.coreg{1}.estimate.ref{1} = ...
    [anatfolders{1} filesep 'mean_functional.img,1'];
jobs{4}.spatial{1}.coreg{1}.estimate.source{1} = anatfiles{1};

% segmentation (normalization estimation)
jobs{4}.spatial{2}.preproc.data{1} = anatfiles{1};

% normalization reslicing
jobs{4}.spatial{3}.normalise{1}.write.subj.matname{1} = ...
    [anatfiles{1}(1:end-6) '_seg_sn.mat'];
jobs{4}.spatial{3}.normalise{1}.write.subj.resample = ...
    cat(1, funcfiles{:}, {[anatfolders{1} filesep 'mean_functional.img,1']});
funcfiles = fnpp(funcfiles, 'w');

% smoothing
jobs{5}.spatial{1}.smooth.data = cat(1, funcfiles{:});
funcfiles = fnpp(funcfiles, 's');

% slice scantiming
if ~opts.st
    jobs(1) = [];
end

% save job?
if ~isempty(opts.jobfile)
end

% run job?
if opts.jobrun
    try
        spm_defaults;
        spm_jobman('run', jobs);
    catch
        rethrow(lasterror);
    end
end

% sub function to prepend variable array with letter
function f = fnpp(f, p)
    if iscell(f)
        for cc = 1:numel(f)
            f{cc} = fnpp(f{cc}, p);
        end
        return;
    end
    [fp{1:3}] = fileparts(f);
    f = [fp{1}, filesep, p, fp{2}, fp{3}];
% end of function fnpp
